##PBAR_FD_Examine.pydatBkgTemp
#
# Make a whole bunch of exploratory plots
#
# 5/1/2013, John Kwong

import matplotlib.pyplot as plt

#####################################
###   GENERAL EXPLORATORY PLOTS   ###
#####################################

# PLOT OF ALL SPECTRA IN SUBPLOT ARRAY

# PLOT ALL SPECTRA FOR A LIST OF DATASETS
groupNameList  = np.array(('ActiveBackground', 'AirDU', 'WaterDU', 'IronDU'))
#groupNameList  = np.array(('ActiveBackground', 'IronDU'))
#groupNameList  = np.array(('WaterDU', 'IronDU'))

# groupNameList  = ['ActiveBackground']

logPlot = 1
filterSpectra = False
normalizeAmplitude = False
normalizeCount = False
subtractBackground = False
subplotx = 4
subploty = 4
subSetIndexList = list(np.arange(8))
subSetIndexList = list(np.arange(5,10))
subSetIndexList = [4]

backgroundIndex = datasetGroupsIndices['ActiveBackground'][3]


# cycle through detectors
for detectorNo in np.arange(40):
    # Generate a new figure window if filled.
    if ((detectorNo % (subplotx * subploty)) == 0):
        f, ax = plt.subplots(subploty, subplotx, sharex='col', sharey='row')
    # Calculate subplot window index
    ii = detectorNo/subploty  - 4*(detectorNo / (subplotx * subploty))
    jj = detectorNo%subplotx
    
    # Cycle through the dataset lists
    for nn in np.arange(len(groupNameList)):
        for dd in subSetIndexList:
            index = datasetGroupsIndices[groupNameList[nn]][dd]  # plot only one
            temp = dat[index][detectorNo,:]

            if subtractBackground:
                temp = temp - dat[backgroundIndex][detectorNo,:]
            if filterSpectra:
                temp = ndimage.filters.gaussian_filter(temp, 1, order = 0)
            if normalizeAmplitude:
                temp = temp / temp[0:50].max()
            if normalizeCount:
                temp = temp / temp[0:250].sum()                        

            ax[ii][jj].plot(temp, label = prefixList[index])
        ax[ii][jj].set_title(str(detectorNo+1))
        
        if ( ii == 0 and jj == 0):
            legend_ = ax[ii][jj].legend()
            # change font size
        if logPlot:
            ax[ii][jj].set_yscale('log')
        else:
            if (not normalizeAmplitude and  not normalizeCount):
                ax[ii][jj].set_ylim((0,20000))            
            ax[ii][jj].set_xlim((0,50))


# SUM SPECTRA

groupNameList  = np.array(('ActiveBackground', 'AirDU', 'WaterDU', 'IronDU'))
groupNameList  = np.array(('ActiveBackground', 'IronDU'))
groupNameList  = np.array(('ActiveBackground', 'AirDU', 'WaterDU', 'IronDU'))

# groupNameList  = ['ActiveBackground']

logPlot = 1
filterSpectra = False
normalizeAmplitude = False
normalizeCount = False
subtractBackground = False

subSetIndexList = list(np.arange(8))
subSetIndexList = list(np.arange(5,10))
subSetIndexList = [1]

backgroundIndex = datasetGroupsIndices['ActiveBackground'][3]
detectorIndices = detGroups['liquid']
#detectorIndices = detGroups['plastic']

plt.figure()
# Cycle through the dataset lists
for nn in np.arange(len(groupNameList)):
    for dd in subSetIndexList:
        index = datasetGroupsIndices[groupNameList[nn]][dd]  # plot only one
        temp = dat[index][detectorIndices,:].sum(axis = 0)
        if subtractBackground:
            temp = temp - dat[backgroundIndex][detectorIndices,:].sum(axis = 0)
        if filterSpectra:
            temp = ndimage.filters.gaussian_filter(temp, 1, order = 0)
        if normalizeAmplitude:
            temp = temp / temp[0:50].max()
        if normalizeCount:
            temp = temp / temp[0:250].sum()
        plot(temp, label = groupNameList[nn] + ', ' + prefixList[index])
        print temp[130:200].sum()
    legend_ = plt.legend()
    if logPlot:
        plt.yscale('log')
        plt.xlim((0, 255))
    else:
        if (not normalizeAmplitude and  not normalizeCount):
            plt.ylim((0,20000))            
        plt.xlim((0,200))
title('Summed Spectra');


# CUMMULATIVE SUM OF SUMMED SPECTRA
groupNameList  = np.array(('ActiveBackground', 'AirDU', 'WaterDU', 'IronDU'))
groupNameList  = np.array(('ActiveBackground', 'IronDU'))
groupNameList  = np.array(('ActiveBackground', 'AirDU', 'WaterDU', 'IronDU'))

# groupNameList  = ['ActiveBackground']

logPlot = 0
filterSpectra = False
normalizeAmplitude = False
normalizeCount = False
subtractBackground = False

subSetIndexList = list(np.arange(8))
subSetIndexList = list(np.arange(5,10))
subSetIndexList = [7]

backgroundIndex = datasetGroupsIndices['ActiveBackground'][3]
detectorIndices = detGroups['liquid']
detectorIndices = detGroups['plastic']

windowBounds = np.array((45, 255))
windowMask = np.zeros(256)
windowMask[windowBounds[0]:(windowBounds[1]+1)] = 1


plt.figure()
# Cycle through the dataset lists
for nn in np.arange(len(groupNameList)):
    for dd in subSetIndexList:
        index = datasetGroupsIndices[groupNameList[nn]][dd]  # plot only one
        temp = dat[index][detectorIndices,:].sum(axis = 0)
        if subtractBackground:
            temp = temp - dat[backgroundIndex][detectorIndices,:].sum(axis = 0)
        if filterSpectra:
            temp = ndimage.filters.gaussian_filter(temp, 1, order = 0)
        if normalizeAmplitude:
            temp = temp / temp[0:50].max()
        if normalizeCount:
            temp = temp / temp[0:250].sum()
        plot(temp / temp.max(), label = groupNameList[nn] + ', ' + prefixList[index])
        plot(  (temp * windowMask ).cumsum() / (temp * windowMask ).sum(), label = groupNameList[nn] + ', ' + prefixList[index])
        print temp[130:180].sum()
    legend_ = plt.legend()
    if logPlot:
        plt.yscale('log')
        plt.xlim((0, 255))
    else:
        #if (not normalizeAmplitude and  not normalizeCount):
            # plt.ylim((0,20000))            
        plt.xlim((0,256))
title('Cummulative Sum of Spectra');


# PLOT SPECTRA FOR A LIST OF DETECTORS FROM A SINGLE DATASET
groupName =  'NoSampleAir'
# groupName =  'NoSampleFe'
# groupName =  'NoSampleH2O'
##groupName = 'NoSampleBackground'
groupName =  'ActiveBackground'

detectorNumber = 1
plt.figure()

index = datasetGroupsIndices[groupName][0]
detectorList = plastic
detectorList = liquid[0:8]
filterSpectra = 0

for ii in range(len(detectorList)):
    i = detectorList[ii]
    if filterSpectra:
        plt.plot(ndimage.filters.gaussian_filter(dat[index][i,:], 4, order = 0) , label = prefixList[index] + ', ' + str(i+1) + ', filtered')
    else:
        plt.plot(dat[index][i,:], label = prefixList[index] + ', ' + str(i+1))

plt.title('detector' + str(detectorNumber))
plt.legend(loc = 1)
plt.grid()
plt.yscale('log')
plt.show()

# PLOT SPECTRA FROM ONE DATASET FOR MANY DATASETS
groupName =  'NoSampleAir'
# groupName =  'NoSampleFe'
# groupName =  'NoSampleH2O'
##groupName = 'NoSampleBackground'
groupName =  'ActiveBackground'

plt.figure()
plt.grid()
datasetIndices = datasetGroupsIndices[groupName]
detectorNumber = 1
filterSpectra = 0

aluminumBounds = np.array((24, 32))

for ii in range(len(datasetIndices)):
    index = datasetIndices[ii]
    
    count = dat[index][detectorNumber,aluminumBounds[0]:aluminumBounds[1]].sum()
    if filterSpectra:
        plt.plot(ndimage.filters.gaussian_filter(dat[index][detectorNumber,:], 4, order = 0), \
                 label = prefixList[index] + ', ' + str(detectorNumber+1) + ', filtered, Al Counts = ' + str(count), \
                 ls = lineStyles[ii/7], color = plotColors[ii%7])
    else:
        plt.plot(dat[index][detectorNumber,:]/count, \
                 label = prefixList[index] + ', ' + str(detectorNumber+1) + ', Al Counts = ' + str(count), \
                 ls = lineStyles[ii/7], color = plotColors[ii%7])

x = np.empty(2)
x.fill(aluminumBounds[0])
y = np.array((0.001, 1e5))
plt.plot(x, y)

x = np.empty(2)
x.fill(aluminumBounds[1])
y = np.array((0.001, 1e5))
plt.plot(x, y)

plt.yscale('log')
plt.title(groupName)
plt.legend(loc = 1)
plt.xlabel('Bin')
plt.ylabel('Counts')
plt.xlim((0, 150))
plt.show()




# FILTERED WAVEFORM AND SLOPE OF FILTERED WAVEFORM

groupName =  'NoSampleAir'
# groupName =  'NoSampleFe'
# groupName =  'NoSampleH2O'
##groupName = 'NoSampleBackground'
detectorNumber = 35
plt.figure()

index = datasetGroupsIndices[groupName][0]
detectorList = plastic
detectorList = liquid[0:8]

for ii in range(len(detectorList)):
    i = detectorList[ii]
    #plt.plot(dat[index][i,:], label = prefixList[i] + ', ' + str(i+1))
    plt.plot(diff(ndimage.filters.gaussian_filter(dat[index][i,:], 4, order = 0)) , label = prefixList[i] + ', ' + str(i+1) + ', filtered')

plt.title('detector' + str(detectorNumber))
plt.legend(loc = 1)
plt.show()
#plt.yscale('log')


# TRY FITTING LINE TO SEVERAL SECTIONS OF THE TAIL

# fit parameters
groupName =  'NoSampleAir'
index = datasetGroupsIndices[groupName][0]
detectorList = liquid
detectorNo = detectorList[5]
detectorNo = 0

binFitBound = np.array([60, 100])
index = 2

# Repeatedly fit lines to the spectrum with different bin windows
pfit = np.zeros((10,2))
bins = np.arange(dat[index][detectorNo,:].shape[0])
jj = 0
for ii in range(0,20,2):
    temp = binFitBound + ii
    cut = (bins > temp[0]) & ( bins < temp[1])
    pfit[jj,:] = np.polyfit(bins[cut], log(dat[index][detectorNo,cut]), 1)
    jj = jj + 1
plt.figure()
plot(bins, log(dat[index][detectorNo,:]))
for ii in range(pfit.shape[0]):
    plot(bins, np.polyval(pfit[ii,:], bins))
    
    
    

# EXAMINE SPECTRA ACROSS DATASETS

#groupName =  'NoSampleAir'
#groupName =  'NoSampleFe'
#groupName =  'NoSampleH2O'
groupName = 'NoSampleBackground'
groupName = 'ActiveBackground'
detectorNumber = 29
detectorNumber = 30

plt.figure()
indices = datasetGroupsIndices[groupName]

for ii in range(len(indices)):
    i = indices[ii]
    plt.plot(dat[i][detectorNumber,:], label = prefixList[i])
plt.title('detector' + str(detectorNumber))
plt.show()
plt.yscale('log')

# EXAMINE SPECTRA ACROSS DATASETS, MATCH MAX AMPLITUDE
groupName = 'ActiveBackground'
#detGroups['plastic'] = plastic
#detGroups['liquid'] = liquid
#detGroups['good'] = np.hstack((np.arange(21), 24))
#detGroups['notgood'] = np.array((22, 23, 24, 26, 27, 28, 29, 34, 36, 38)) - 1
#detGroups['bad'] = np.array((30, 31, 32, 33, 35, 39)) - 1
#detGroups['broken'] = 36

detectorNo = detGroups['good'][0]
#detectorNo = detGroups['notgood'][0]
#detectorNo = detGroups['bad'][0]

#detectorNo = 30
filterSpectra = True
normalizeAmplitude = False
normalizeCount = False
plt.figure()

indicesList = datasetGroupsIndices[groupName]

for ii in range(len(indicesList)):
    # dataset index
    i = indicesList[ii]
    # unmodified spectrum
    temp = dat[i][detectorNo,:]
    # modify the spectrum
    if filterSpectra:
        temp = ndimage.filters.gaussian_filter(temp, 1, order = 0)
    if normalizeAmplitude:
        temp = temp / temp[0:50].max()
    if normalizeCount:
        temp = temp / temp[0:250].sum()
    plt.plot(temp, label = (prefixList[i] + ', ' + datasetDescription[i]))
plt.title('detector' + str(detectorNo+1))
plt.show()
#plt.yscale('log')
plt.legend(loc = 1)
